Introduction

Background Information

Honeybees are very important organisms to the world at large. Honeybees are important pollinators of flowering plants, and thus are essential to the agriculture industry. Despite their vital importance, honeybee populations have been declining for nearly thirty years. One of the main reasons for this decline in honeybee numbers is the Varroa mite (Varroa destructor), an external parasite of honeybees. Varroa mites are responsible for causing Varroosis in honeybees, weakening adult bees by feeding on their haemolymph and thereby increasing their susceptibility to other diseases. The current protocol in place to help rid honeybee colonies of Varroa mites is through miticides. However, this method of pest control is not sustainable as the use of miticides may further weaken the honeybees or contaminate hive products. More importantly, an increase in resistance to miticides by Varroa mites has been seen. Consequently, it is imperative to develop a new means of controlling these pests to increase honeybee populations in the long term.
Varroa mites on adult male honey bee (Apis melifera) (genetic literacy project.org)

Research Objectives

This research aims to examine various odorants and their effects on both Varroa mites and honeybees to analyze their attractant or repellent properties on the two species. Additionally, the effects of the odorants on honeybee and Varroa mite behaviour, as well as effects on in-hive production, will be examined.

The outcome of this research project will be to implement the use of odorants as a more sustainable means of controlling Varroa mite populations, without the use of toxic pesticides. Successful implementation of the results of this study could help slow the decline of honeybee populations worldwide by decreasing the presence of Varroa mites in honeybee colonies.

Progress

Over the summer, Mike conducted electrotarsograms (ETGs) on Varroa mites to examine their physiological sensitivity to several previously identified attractive odourants.

ETGs on Varroa mites were conducted using single odourants previously identified as evoking a response in Varroa mites, and were tested in four different concentrations (10ng, 1ng, 100μg, 10μg). These are typical concentrations used in ETG studies and often encompass concentrations close to what Varroa mites would encounter in their natural environment. Relative to a control stimulus, these reactions can be used to infer the relative importance of each odourant to Varroa mites in the hive environment.

Questions

The questions we want to answer from the dataset are as follows:
1. Are there mite trials that are outliers?
2. Are there odours which consistantly evoke a response?
3. Are there concentration-dependent trends?
4. Can we quantify and account for between-trial (between-mite) variability?

Loading Data

my.df <- read_csv('DataForR.csv')
## Parsed with column specification:
## cols(
##   trial = col_integer(),
##   odour = col_character(),
##   conc = col_double(),
##   percent = col_double(),
##   gb = col_character(),
##   normalized = col_double()
## )

Look at Data

str(my.df)
## Classes 'tbl_df', 'tbl' and 'data.frame':    1165 obs. of  6 variables:
##  $ trial     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ odour     : chr  "2-Hydroxyhexanoic Acid (DCM)" "2-Hydroxyhexanoic Acid (DCM)" "2-Hydroxyhexanoic Acid (DCM)" "Methyl Palmitate" ...
##  $ conc      : num  0.01 0.1 1 0.01 0 10 0.01 1 0.01 0.1 ...
##  $ percent   : num  53.6 142.5 157.9 137.7 100 ...
##  $ gb        : chr  "Good" "Good" "Good" "Good" ...
##  $ normalized: num  13.3 129.8 136.7 127.4 100 ...
##  - attr(*, "spec")=List of 2
##   ..$ cols   :List of 6
##   .. ..$ trial     : list()
##   .. .. ..- attr(*, "class")= chr  "collector_integer" "collector"
##   .. ..$ odour     : list()
##   .. .. ..- attr(*, "class")= chr  "collector_character" "collector"
##   .. ..$ conc      : list()
##   .. .. ..- attr(*, "class")= chr  "collector_double" "collector"
##   .. ..$ percent   : list()
##   .. .. ..- attr(*, "class")= chr  "collector_double" "collector"
##   .. ..$ gb        : list()
##   .. .. ..- attr(*, "class")= chr  "collector_character" "collector"
##   .. ..$ normalized: list()
##   .. .. ..- attr(*, "class")= chr  "collector_double" "collector"
##   ..$ default: list()
##   .. ..- attr(*, "class")= chr  "collector_guess" "collector"
##   ..- attr(*, "class")= chr "col_spec"
summary(my.df)
##      trial          odour                conc           percent       
##  Min.   : 1.00   Length:1165        Min.   : 0.000   Min.   :  17.59  
##  1st Qu.: 9.00   Class :character   1st Qu.: 0.010   1st Qu.:  78.07  
##  Median :13.00   Mode  :character   Median : 0.100   Median : 100.00  
##  Mean   :12.22                      Mean   : 2.665   Mean   : 117.17  
##  3rd Qu.:16.00                      3rd Qu.: 1.000   3rd Qu.: 125.02  
##  Max.   :18.00                      Max.   :10.000   Max.   :2507.58  
##       gb              normalized     
##  Length:1165        Min.   :-368.60  
##  Class :character   1st Qu.:  71.91  
##  Mode  :character   Median : 100.00  
##                     Mean   :  89.80  
##                     3rd Qu.: 120.01  
##                     Max.   : 196.01
hist(my.df$normalized, 
     xlab="Normalized ETG Response", ylab="Frequency", 
     main="Frequency of Normalized ETG Responses")

p <- ggplot(data=my.df, aes(normalized))
p + geom_density()

Doing this produces a histogram of the percent reaction by Varroa mites. This is an easy way to see how the data is distributed.

Changing Trial number and Concentration to Factors

#change integer to factor for trials
my.df$trial<-factor(my.df$trial)

#lists the levels of factors
levels(my.df$trial)
##  [1] "1"  "2"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14" "15"
## [15] "16" "17" "18"
#change number to factor for concentration "conc"
my.df$conc<-factor(my.df$conc)
levels(my.df$conc)
## [1] "0"    "0.01" "0.1"  "1"    "10"
# types of odours
my.df$odour<-factor(my.df$odour)
levels(my.df$odour)
##  [1] "1-Hexadecanol"                "2-Heptanol"                  
##  [3] "2-Heptanone"                  "2-Hydroxyhexanoic Acid (DCM)"
##  [5] "2-Nonanal"                    "Benzoic Acid"                
##  [7] "Butyric Acid"                 "Dodecyl Aldehyde"            
##  [9] "Ethyl Palmitate"              "Geraniol"                    
## [11] "Heptacosane"                  "Heptadecane"                 
## [13] "Hexadecanal"                  "Hexane"                      
## [15] "Methyl Oleate"                "Methyl Palmitate"            
## [17] "Octadecanol"                  "Octanoic Acid"               
## [19] "Palmitic Acid"                "trans-Nerolidol"

Creating a Faceted Plot for All Odours at All Concentrations

Step1: Filter Out Bad Runs

We discovered that these runs were highly variable and produced messy plots, so we will remove them for now.

good_runs<- filter(my.df,
                    !(odour %in% c("1-Hexadecanol", "Methyl Palmitate", "Hexane",   "Palmitic Acid" )))

Step2: Normalize Scores at Zero Instead of 100

subsetted <- mutate(good_runs,
                    Zresponse = normalized - 100)

Look at Normalized Data Distribution

hist(subsetted$Zresponse,xlab="Normalized ETG Response",ylab="Frequency", main="Frequency of Normalized ETG Responses")

Step3: Filter Out Negative Responses

A negative response makes no sense physiologically, so this is an indication of a bad run. We remove these negative values from analysis.

nonegatives <- subsetted %>%
  mutate(Zresponse= ifelse(Zresponse <0, NA, Zresponse)) 
str(nonegatives)
## Classes 'tbl_df', 'tbl' and 'data.frame':    706 obs. of  7 variables:
##  $ trial     : Factor w/ 17 levels "1","2","4","5",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ odour     : Factor w/ 20 levels "1-Hexadecanol",..: 4 4 4 15 12 12 12 2 5 5 ...
##  $ conc      : Factor w/ 5 levels "0","0.01","0.1",..: 2 3 4 5 2 3 4 3 2 3 ...
##  $ percent   : num  53.6 142.5 157.9 202.7 44.1 ...
##  $ gb        : chr  "Good" "Good" "Good" "Good" ...
##  $ normalized: num  13.3 129.8 136.7 150.7 -27 ...
##  $ Zresponse : num  NA 29.8 36.7 50.7 NA ...

Step4: Make the Plot

Now that we have manipulated the data in a way that is useful to us, we will make our faceted plot.

ggplot(nonegatives)+
  geom_boxplot(aes(x=conc,y= Zresponse, group= conc))+
  xlab("Concentration")+ylab("Normalized Percent Reaction")+ggtitle("Faceted ETG Plots")+
  theme_bw(10)+
  facet_wrap(~odour)
## Warning: Removed 366 rows containing non-finite values (stat_boxplot).

From this plot, we can see the overall trend in ETG responses from Varroa mites at each concentration for each compound tested.

Step5: Taking temperature and humidity into account

We are still waiting on temperature and humidity data, but Mike has been in touch with someone regarding this.

absoluteweather<- read.csv('absoluteweather.csv')

library(tidyverse)
absoluteweather$conc <- factor(absoluteweather$conc)
absoluteweather$trial <- factor(absoluteweather$trial)
str(absoluteweather)
## 'data.frame':    1165 obs. of  13 variables:
##  $ trial     : Factor w/ 17 levels "1","2","4","5",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ odour     : Factor w/ 20 levels "1-Hexadecanol",..: 4 4 4 16 14 15 1 16 12 12 ...
##  $ conc      : Factor w/ 5 levels "0","0.01","0.1",..: 2 3 4 2 1 5 2 4 2 3 ...
##  $ absamp    : num  3.02 8.09 9 7.95 5.91 ...
##  $ linearhex : num  5.64 5.67 5.7 5.77 5.91 ...
##  $ percent   : num  53.6 142.5 157.9 137.7 100 ...
##  $ normalized: num  13.3 129.8 136.7 127.4 100 ...
##  $ gb        : Factor w/ 2 levels "Bad","Good": 2 2 2 2 2 2 2 2 2 2 ...
##  $ humid     : num  69 69 69.1 69.2 69.3 ...
##  $ temp      : num  20.8 20.8 20.8 20.7 20.7 ...
##  $ dewp      : num  14.9 14.9 14.9 14.9 14.9 ...
##  $ hour      : num  2013 2014 2014 2016 2019 ...
##  $ min       : num  13.1 13.8 14.5 16 18.9 ...
ggplot(data = absoluteweather, mapping = aes(humid, absamp)) + 
  geom_jitter(aes(colour = absoluteweather$trial), width = 0.25) +
  ggtitle("Amplitude vs Humidity")

ggplot(data = absoluteweather, mapping = aes(temp, absamp)) + 
  geom_jitter(aes(colour = absoluteweather$trial), width = 0.25) +
  ggtitle("Amplitude vs Temperature")

ggplot(data = absoluteweather, mapping = aes(dewp, absamp)) + 
  geom_jitter(aes(colour = absoluteweather$trial), width = 0.25) +
  ggtitle("Amplitude vs Dewpoint")

ggplot(data = absoluteweather, mapping = aes(absamp))+
  geom_histogram(bins = 50)

# correlations
cor(absoluteweather$absamp, absoluteweather$temp)
## [1] -0.4868699
cov(absoluteweather$absamp, absoluteweather$temp)
## [1] -234.6016
corr_amp_tempr <- cor.test(x=absoluteweather$absamp, 
                           y=absoluteweather$temp, method = 'spearman')
## Warning in cor.test.default(x = absoluteweather$absamp, y = absoluteweather
## $temp, : Cannot compute exact p-value with ties
corr_amp_tempr
## 
##  Spearman's rank correlation rho
## 
## data:  absoluteweather$absamp and absoluteweather$temp
## S = 367540000, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.3946761
corr_amp_humid <- cor.test(x=absoluteweather$absamp, 
                           y=absoluteweather$humid, method = 'spearman')
## Warning in cor.test.default(x = absoluteweather$absamp, y = absoluteweather
## $humid, : Cannot compute exact p-value with ties
corr_amp_humid
## 
##  Spearman's rank correlation rho
## 
## data:  absoluteweather$absamp and absoluteweather$humid
## S = 103650000, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.606672
corr_amp_dewp<- cor.test(x=absoluteweather$absamp, 
                         y=absoluteweather$dewp, method = 'spearman')
## Warning in cor.test.default(x = absoluteweather$absamp, y = absoluteweather
## $dewp, : Cannot compute exact p-value with ties
corr_amp_dewp
## 
##  Spearman's rank correlation rho
## 
## data:  absoluteweather$absamp and absoluteweather$dewp
## S = 195510000, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.2580989
library(fifer)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
spearman.plot(absoluteweather$absamp, absoluteweather$temp)

weather_matrix <- data.matrix(absoluteweather[7:10], rownames.force = NA)
pairs(weather_matrix) # plots all the stuff simultaneiously

Plotting absolute amplitude (corrected to Hexane blank) versus time of trial.

absoluteweather <- read.csv('absoluteweather.csv')
str(absoluteweather)
## 'data.frame':    1165 obs. of  13 variables:
##  $ trial     : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ odour     : Factor w/ 20 levels "1-Hexadecanol",..: 4 4 4 16 14 15 1 16 12 12 ...
##  $ conc      : num  0.01 0.1 1 0.01 0 10 0.01 1 0.01 0.1 ...
##  $ absamp    : num  3.02 8.09 9 7.95 5.91 ...
##  $ linearhex : num  5.64 5.67 5.7 5.77 5.91 ...
##  $ percent   : num  53.6 142.5 157.9 137.7 100 ...
##  $ normalized: num  13.3 129.8 136.7 127.4 100 ...
##  $ gb        : Factor w/ 2 levels "Bad","Good": 2 2 2 2 2 2 2 2 2 2 ...
##  $ humid     : num  69 69 69.1 69.2 69.3 ...
##  $ temp      : num  20.8 20.8 20.8 20.7 20.7 ...
##  $ dewp      : num  14.9 14.9 14.9 14.9 14.9 ...
##  $ hour      : num  2013 2014 2014 2016 2019 ...
##  $ min       : num  13.1 13.8 14.5 16 18.9 ...
absoluteweather$trial <- as.factor(absoluteweather$trial)
absoluteweather$conc <- as.factor(absoluteweather$conc)

Step1: Normalize scores at zero

subsetted <- mutate(absoluteweather,
                    Zresponse = normalized - 100)

Step2: Filter out negative responses

#Not sure whether to replace nonsense scores with NA or 0
nonegatives <- subsetted %>%
  mutate(Zresponse= ifelse(Zresponse <0, NA, Zresponse)) 
str(nonegatives$Zresponse)
##  num [1:1165] NA 29.8 36.7 27.4 0 ...

ETG amplitude plots faceted by trial

ggplot(data = nonegatives, aes(min, absamp))+
  geom_jitter(aes(colour = nonegatives$odour), width = 0.25)+
  xlab("Time")+ylab("Absamp")+ggtitle("Faceted ETG Plots")+
  theme_bw(10)+
  facet_wrap(~trial)

Cut out the Noisy Trails

nonegatives$trial <- as.numeric(nonegatives$trial)
goodtrials <- nonegatives[nonegatives$trial %in% c(1:8,11, 12), ]
goodtrials$trial <- as.factor(goodtrials$trial)
str(goodtrials)
## 'data.frame':    530 obs. of  14 variables:
##  $ trial     : Factor w/ 10 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ odour     : Factor w/ 20 levels "1-Hexadecanol",..: 4 4 4 16 14 15 1 16 12 12 ...
##  $ conc      : Factor w/ 5 levels "0","0.01","0.1",..: 2 3 4 2 1 5 2 4 2 3 ...
##  $ absamp    : num  3.02 8.09 9 7.95 5.91 ...
##  $ linearhex : num  5.64 5.67 5.7 5.77 5.91 ...
##  $ percent   : num  53.6 142.5 157.9 137.7 100 ...
##  $ normalized: num  13.3 129.8 136.7 127.4 100 ...
##  $ gb        : Factor w/ 2 levels "Bad","Good": 2 2 2 2 2 2 2 2 2 2 ...
##  $ humid     : num  69 69 69.1 69.2 69.3 ...
##  $ temp      : num  20.8 20.8 20.8 20.7 20.7 ...
##  $ dewp      : num  14.9 14.9 14.9 14.9 14.9 ...
##  $ hour      : num  2013 2014 2014 2016 2019 ...
##  $ min       : num  13.1 13.8 14.5 16 18.9 ...
##  $ Zresponse : num  NA 29.8 36.7 27.4 0 ...
remove(subsetted)

Histogram of the normalized data

ggplot(data = nonegatives, mapping = aes(Zresponse))+
  geom_histogram(bins = 50)
## Warning: Removed 577 rows containing non-finite values (stat_bin).

Analysis of Variance (ANOVA)

plot(aov_out <- aov(Zresponse~odour*conc*trial, data = goodtrials))
## Warning: not plotting observations with leverage one:
##   1, 2, 5, 8, 12, 13, 14, 19, 21, 22, 23, 24, 26, 28, 31, 32, 33, 38, 39, 40, 43, 45, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 64, 65, 67, 68, 69, 70, 71, 85, 86, 87, 88, 89, 90, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 116, 117, 118, 119, 120, 125, 126, 127, 128, 129, 131, 132, 133, 136, 137, 140, 141, 146, 154, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 171, 172, 177, 178, 179, 180, 181, 183, 184, 185, 186, 187, 191, 192, 193, 194, 195, 196, 197, 198, 199, 205, 206, 207, 208, 209, 210, 211, 212, 217, 218, 219, 222, 230, 233, 238, 239, 240, 241, 242, 243, 244, 251, 252, 253, 254, 257, 258, 260, 267, 268, 269

## Warning: not plotting observations with leverage one:
##   1, 2, 5, 8, 12, 13, 14, 19, 21, 22, 23, 24, 26, 28, 31, 32, 33, 38, 39, 40, 43, 45, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 64, 65, 67, 68, 69, 70, 71, 85, 86, 87, 88, 89, 90, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 116, 117, 118, 119, 120, 125, 126, 127, 128, 129, 131, 132, 133, 136, 137, 140, 141, 146, 154, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 171, 172, 177, 178, 179, 180, 181, 183, 184, 185, 186, 187, 191, 192, 193, 194, 195, 196, 197, 198, 199, 205, 206, 207, 208, 209, 210, 211, 212, 217, 218, 219, 222, 230, 233, 238, 239, 240, 241, 242, 243, 244, 251, 252, 253, 254, 257, 258, 260, 267, 268, 269

## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced

summary(aov_out)
##                  Df Sum Sq Mean Sq F value   Pr(>F)    
## odour            19  39802    2095   6.955 8.29e-11 ***
## conc              3    523     174   0.579   0.6305    
## trial             9  40810    4534  15.054 1.64e-14 ***
## odour:conc       41  11033     269   0.893   0.6497    
## odour:trial      55  23456     426   1.416   0.0722 .  
## conc:trial       25   6152     246   0.817   0.7108    
## odour:conc:trial 31   7147     231   0.765   0.7978    
## Residuals        88  26507     301                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 258 observations deleted due to missingness

Messing Around With Models

Linear Model (LM)

library(tidyverse)
library(modelr)
options(na.action = na.warn)
ggplot(goodtrials, aes(trial, absamp))+
  geom_point(aes(colour=odour))

Zresponse ~ concentration+odour

mod01 <- lm(Zresponse~conc+odour, data = goodtrials)
## Warning: Dropping 258 rows with missing values
#fits red dots as predictions according to the model
grid <- goodtrials %>% 
  data_grid(odour,conc) %>% 
  add_predictions(mod01)
## Warning in predict.lm(model, data): prediction from a rank-deficient fit
## may be misleading
grid
## # A tibble: 100 x 3
##    odour         conc                 pred
##    <fct>         <fct>               <dbl>
##  1 1-Hexadecanol 0       0.000000000000165
##  2 1-Hexadecanol 0.01   27.2              
##  3 1-Hexadecanol 0.1    29.9              
##  4 1-Hexadecanol 1      30.9              
##  5 1-Hexadecanol 10     27.6              
##  6 2-Heptanol    0     - 7.80             
##  7 2-Heptanol    0.01   19.4              
##  8 2-Heptanol    0.1    22.1              
##  9 2-Heptanol    1      23.1              
## 10 2-Heptanol    10     19.8              
## # ... with 90 more rows
#plots our data
ggplot(goodtrials, aes(x = conc)) + 
  geom_point(aes(y = Zresponse)) +
  geom_point(data = grid, aes(y = pred), colour = "red", size = 0.2)+
  facet_wrap (~odour)
## Warning: Removed 258 rows containing missing values (geom_point).

Assessment of the Linear Model

This model doesn’t appear to fit the data very well

Poisson GLM

#I tells the below function to change Zresponse to an integer 
mod01 <- glm(I(as.integer(goodtrials$Zresponse))~conc+odour, data = goodtrials, poisson)
## Warning: Dropping 258 rows with missing values
## Warning: glm.fit: algorithm did not converge
par(mfrow=c(2, 2))
summary(mod01)
## 
## Call:
## glm(formula = I(as.integer(goodtrials$Zresponse)) ~ conc + odour, 
##     family = poisson, data = goodtrials)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -41.048   -5.423   -1.785    1.715    9.571  
## 
## Coefficients: (1 not defined because of singularities)
##                                    Estimate Std. Error  z value Pr(>|z|)
## (Intercept)                        6.736328   0.003304 2038.804  < 2e-16
## conc0.01                          -3.454085   0.095869  -36.029  < 2e-16
## conc0.1                           -3.363648   0.100125  -33.594  < 2e-16
## conc1                             -3.332014   0.098638  -33.780  < 2e-16
## conc10                            -3.442636   0.100819  -34.147  < 2e-16
## odour2-Heptanol                   -0.307637   0.118917   -2.587 0.009682
## odour2-Heptanone                   0.234648   0.107940    2.174 0.029715
## odour2-Hydroxyhexanoic Acid (DCM) -0.080565   0.111600   -0.722 0.470348
## odour2-Nonanal                     0.024033   0.109300    0.220 0.825966
## odourBenzoic Acid                 -0.608394   0.123951   -4.908 9.18e-07
## odourButyric Acid                 -0.642706   0.178872   -3.593 0.000327
## odourDodecyl Aldehyde              0.050435   0.127273    0.396 0.691904
## odourEthyl Palmitate               0.265458   0.138084    1.922 0.054550
## odourGeraniol                     -0.785819   0.148409   -5.295 1.19e-07
## odourHeptacosane                  -0.130434   0.126651   -1.030 0.303071
## odourHeptadecane                   0.185362   0.106257    1.744 0.081079
## odourHexadecanal                   0.278775   0.110403    2.525 0.011568
## odourHexane                              NA         NA       NA       NA
## odourMethyl Oleate                 0.274981   0.106797    2.575 0.010030
## odourMethyl Palmitate              0.116862   0.098339    1.188 0.234691
## odourOctadecanol                  -0.198848   0.157146   -1.265 0.205739
## odourOctanoic Acid                 0.034727   0.110097    0.315 0.752441
## odourPalmitic Acid                 0.598129   0.174819    3.421 0.000623
## odourtrans-Nerolidol              -0.231625   0.177177   -1.307 0.191108
##                                      
## (Intercept)                       ***
## conc0.01                          ***
## conc0.1                           ***
## conc1                             ***
## conc10                            ***
## odour2-Heptanol                   ** 
## odour2-Heptanone                  *  
## odour2-Hydroxyhexanoic Acid (DCM)    
## odour2-Nonanal                       
## odourBenzoic Acid                 ***
## odourButyric Acid                 ***
## odourDodecyl Aldehyde                
## odourEthyl Palmitate              .  
## odourGeraniol                     ***
## odourHeptacosane                     
## odourHeptadecane                  .  
## odourHexadecanal                  *  
## odourHexane                          
## odourMethyl Oleate                *  
## odourMethyl Palmitate                
## odourOctadecanol                     
## odourOctanoic Acid                   
## odourPalmitic Acid                ***
## odourtrans-Nerolidol                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance:  6453.5  on 271  degrees of freedom
## Residual deviance: 71316.0  on 249  degrees of freedom
##   (258 observations deleted due to missingness)
## AIC: 72458
## 
## Number of Fisher Scoring iterations: 25
plot(mod01)
## Warning: not plotting observations with leverage one:
##   520
## Warning: not plotting observations with leverage one:
##   520

Assessment of the Poisson GLM

This model doesn’t appear to fit the data well

Gaussian GLM

mod02 <- glm(goodtrials$Zresponse~conc+odour, data = goodtrials)
## Warning: Dropping 258 rows with missing values
par(mfrow=c(2, 2))
summary(mod02)
## 
## Call:
## glm(formula = goodtrials$Zresponse ~ conc + odour, data = goodtrials)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -36.056  -14.001   -0.068    9.043   65.476  
## 
## Coefficients: (1 not defined because of singularities)
##                                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)                        1.654e-13  3.400e+00   0.000  1.00000
## conc0.01                           2.724e+01  1.132e+01   2.406  0.01684
## conc0.1                            2.993e+01  1.184e+01   2.528  0.01209
## conc1                              3.091e+01  1.168e+01   2.647  0.00864
## conc10                             2.757e+01  1.189e+01   2.320  0.02118
## odour2-Heptanol                   -7.799e+00  1.290e+01  -0.604  0.54613
## odour2-Heptanone                   7.476e+00  1.254e+01   0.596  0.55152
## odour2-Hydroxyhexanoic Acid (DCM) -2.186e+00  1.250e+01  -0.175  0.86132
## odour2-Nonanal                     5.983e-01  1.240e+01   0.048  0.96154
## odourBenzoic Acid                 -1.290e+01  1.273e+01  -1.014  0.31172
## odourButyric Acid                 -1.326e+01  1.650e+01  -0.804  0.42219
## odourDodecyl Aldehyde              1.108e+00  1.457e+01   0.076  0.93942
## odourEthyl Palmitate               7.835e+00  1.660e+01   0.472  0.63739
## odourGeraniol                     -1.615e+01  1.410e+01  -1.145  0.25323
## odourHeptacosane                  -3.669e+00  1.397e+01  -0.263  0.79311
## odourHeptadecane                   5.537e+00  1.223e+01   0.453  0.65110
## odourHexadecanal                   8.822e+00  1.297e+01   0.680  0.49700
## odourHexane                               NA         NA      NA       NA
## odourMethyl Oleate                 8.774e+00  1.243e+01   0.706  0.48102
## odourMethyl Palmitate              3.294e+00  1.113e+01   0.296  0.76752
## odourOctadecanol                  -5.052e+00  1.660e+01  -0.304  0.76114
## odourOctanoic Acid                 7.602e-01  1.254e+01   0.061  0.95169
## odourPalmitic Acid                 2.171e+01  2.433e+01   0.892  0.37322
## odourtrans-Nerolidol              -6.257e+00  1.865e+01  -0.336  0.73750
##                                     
## (Intercept)                         
## conc0.01                          * 
## conc0.1                           * 
## conc1                             **
## conc10                            * 
## odour2-Heptanol                     
## odour2-Heptanone                    
## odour2-Hydroxyhexanoic Acid (DCM)   
## odour2-Nonanal                      
## odourBenzoic Acid                   
## odourButyric Acid                   
## odourDodecyl Aldehyde               
## odourEthyl Palmitate                
## odourGeraniol                       
## odourHeptacosane                    
## odourHeptadecane                    
## odourHexadecanal                    
## odourHexane                         
## odourMethyl Oleate                  
## odourMethyl Palmitate               
## odourOctadecanol                    
## odourOctanoic Acid                  
## odourPalmitic Acid                  
## odourtrans-Nerolidol                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 462.2682)
## 
##     Null deviance: 155429  on 271  degrees of freedom
## Residual deviance: 115105  on 249  degrees of freedom
##   (258 observations deleted due to missingness)
## AIC: 2464.9
## 
## Number of Fisher Scoring iterations: 2
plot(mod02)
## Warning: not plotting observations with leverage one:
##   520
## Warning: not plotting observations with leverage one:
##   520

Assessment of the Gaussian GLM

This model appears to fit the data better than the poisson model

Binomial GLM

bindf <- goodtrials %>%
  mutate(BinZresponse= ifelse(Zresponse > 0, 1, 0))
mod03 <- glm(bindf$BinZresponse~conc+odour, data = bindf, binomial)
## Warning: Dropping 258 rows with missing values
## Warning: glm.fit: algorithm did not converge
par(mfrow=c(2, 2))
summary(mod03)
## 
## Call:
## glm(formula = bindf$BinZresponse ~ conc + odour, family = binomial, 
##     data = bindf)
## 
## Deviance Residuals: 
##        Min          1Q      Median          3Q         Max  
## -2.409e-06   2.409e-06   2.409e-06   2.409e-06   2.409e-06  
## 
## Coefficients: (1 not defined because of singularities)
##                                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)                       -2.657e+01  5.631e+04       0        1
## conc0.01                           5.313e+01  1.875e+05       0        1
## conc0.1                            5.313e+01  1.961e+05       0        1
## conc1                              5.313e+01  1.935e+05       0        1
## conc10                             5.313e+01  1.969e+05       0        1
## odour2-Heptanol                    2.117e-08  2.137e+05       0        1
## odour2-Heptanone                   2.471e-09  2.077e+05       0        1
## odour2-Hydroxyhexanoic Acid (DCM)  6.009e-09  2.071e+05       0        1
## odour2-Nonanal                     2.924e-09  2.053e+05       0        1
## odourBenzoic Acid                 -2.023e-09  2.109e+05       0        1
## odourButyric Acid                 -4.927e-06  2.732e+05       0        1
## odourDodecyl Aldehyde              9.531e-10  2.413e+05       0        1
## odourEthyl Palmitate              -3.604e-08  2.750e+05       0        1
## odourGeraniol                      3.838e-08  2.336e+05       0        1
## odourHeptacosane                  -9.987e-09  2.314e+05       0        1
## odourHeptadecane                  -1.393e-09  2.025e+05       0        1
## odourHexadecanal                  -4.719e-10  2.148e+05       0        1
## odourHexane                               NA         NA      NA       NA
## odourMethyl Oleate                 1.646e-09  2.059e+05       0        1
## odourMethyl Palmitate              2.710e-09  1.844e+05       0        1
## odourOctadecanol                  -3.604e-08  2.750e+05       0        1
## odourOctanoic Acid                -4.931e-06  2.076e+05       0        1
## odourPalmitic Acid                -1.437e-08  4.030e+05       0        1
## odourtrans-Nerolidol               2.765e-08  3.089e+05       0        1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2.2716e+02  on 271  degrees of freedom
## Residual deviance: 1.5780e-09  on 249  degrees of freedom
##   (258 observations deleted due to missingness)
## AIC: 46
## 
## Number of Fisher Scoring iterations: 25
plot(mod03)
## Warning: not plotting observations with leverage one:
##   520
## Warning: not plotting observations with leverage one:
##   520

Assessment of the Binomial GLM

This model appears to fit the data terribly

Back to top